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f^*) ■ Abstract 

It is shown that F(R)-modified gravitational theories lead to curvature oscillations in 
astrophysical systems with rising energy density. The frequency and the amplitude of such 
oscillations could be very high and would lead to noticeable production of energetic cosmic 
ray particles. 



There are two popular phenomenological explanations of the observed accelerated universe 
expansion. It is either prescribed to existence of the so called dark energy (DE) with the equation 



of state: 



P*-Q, (1) 



where P and g are respectively pressure and energy densities. Another way to describe cosmo- 
logical acceleration is a modification of the classical action of the general relativity (GR) by an 
addition of a non-linear function of curvature, F(R): 
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where function F(R) is chosen in such a way that the modified GR equations have a solution 
R = const in absence of any matter source. Taking the trace g^,uSA/5g^ u = 0, we find: 

W i F > R _ R + RF ' R _2F = 8nTH/m 2 Pl . (3) 

where T> 2 = V^D^ is the covariant D'Alambertian operator, F' R = dF/dR, and T^ u is the energy- 
momentum tensor of matter. 
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To describe the astronomical data about cosmological acceleration the solution of the equation 



R c - R C F' R {R C ) + 2F(R C ) = (4) 

should be equal to R c = —?>2-kQ<\q c / m 2 pl , where £l\ ps 0.75 is the fraction of the observed vacuum- 
like energy density and g c ~ 10 -29 g/cm 3 is the total cosmological energy density. 

The pioneering suggestion [1] of gravity modification with F(R) ~ fJ. 4 /R suffered from strong 
instabilities in celestial bodies [2]. To cure these instability further modifications of gravity have 
been suggested [3H5]- These and some other versions are reviewed e.g. in refs. [6j[7]. The 
suggested modifications, however, may lead to infinite-i? singularities in the past cosmological 
history [6] and in the future in astronomical systems with rising energy/matter density [HJE]. 
These singularities can be successfully eliminated by an addition of the i2 2 -term into the action. 
Such a contribution naturally appears as a result of quantum corrections due to matter loops 
in curved space-time [lOUllj . Possibly in astrophysical systems singularity may be avoided even 
without R 2 if one takes into account distortion of the background flat Minkowsky space-time by 
an impact of rising R. However, it would take place if R very much differs from its GR value and 
such a deviation from normal gravity would be observable. 

In what follows we consider the version of the modified gravity suggested by Starobinsky [5]: 



F(R) = -XR 



(5) 



6m 2 



where n is an integer, A > 0, and |i?o| ~ 1/^7 > wnere tu ~ 13 Gyr is the universe age. Parameter 
m is bounded from below, m > 10 5 GeV, to preserve successful predictions of BBN, see e.g. [14] . 

Other models j3|4] demonstrate similar behavior. We study the evolution of R in a contracting 
astrophysical system and show that R oscillates around its GR value with quite a large ampli- 
tude, efficiently producing elementary particles with masses below the oscillation frequency. The 
particle production damps the oscillations and may suppress or even eliminate the singularity, 
which would appear if R 2 term was absent. 

Cosmology with R 2 term and in particular the particle production by the oscillating curvature 
was studied in the early works |11H13| . There was renewed interest to this problem [14^115] 
stimulated by the study of the interplay of the infrared and ultraviolet terms in eq. ([5]) . 

In this paper we study solutions of eq. (|3|) in the system of non-relativistic particles with 
rising energy density which we approximate as g = Qo(l + t/t contr ), where t con t r is the effective 
time of contraction. We assume that the gravity of matter is not strong and thus the background 
metric can be considered as flat. Written in terms of R equation of motion ([3]) contains non- linear 
derivative terms: 

[(d t R) 2 - (djRfj/R, (6) 

which makes it very inconvenient for qualitative analysis. To avoid that, we introduce, instead 
of R, the new function, proportional to F'(R): 



2 An \RoJ 

where To = Sir qq / m 2 pl . It is also convenient to introduce dimensionless quantities z = g(t)/go, 
y = —R/Tq, g = l/[Q\n(mtjj) 2 ] (g m o / g c ) 2n+2 , and dimensionless time r = mt^Jg. In terms of 
these quantities and in the limit of R 3> -Ro and sufficiently homogeneous ^ (both conditions are 
naturally fulfilled) the equation of motion for £ takes a very simple form: 

£" + dU/d£ = , (8) 
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where prime denotes derivative with respect to r and dll /d£ = z — ?/(£). Unfortunately y cannot 
be expressed through £ analytically. We have only inverse relation 

£ = l/y 2n+1 - gy, (9) 

so we have to use different approximate expressions in different ranges of £. 

In what follows we make analytical estimates of the amplitude of the curvature oscillations 
and compare them with numerical calculations for some values of the parameters. The agreement 
is generally quite good. However, we cannot do numerical calculations for realistically tiny values 
of g or huge frequencies of oscillations due to numerical instability and/or too long computational 
time. So an agreement of the analytical results with numerical ones where the later can be done 
allows to conclude that analytical estimates can be valid for high frequency and small g as well. 

The minimum of potential U(£) is located at = z(t), so it moves with time according to 

Uin(r) = z(T)-^ n+ ^-gz(T). (10) 

It is intuitively clear that even if initially £ takes its GR value £ = £, m in it would not catch the 
motion of the minimum and as a result it starts to oscillate around it. Dimensionless frequency 
of small oscillations, Q, is determined by: 



y=z 



2n + l 

z 2n+2 
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Note that physical frequency is oj = CI m^Jg. 

For small oscillations £ = £ m i„, + £i takes the form: 

£(r) « Uin(r) + a(r) sin[F(r) + 6} (12) 

where 5 is a phase determined by the initial conditions and 

F(t)= f dT'fl(r'). (13) 



For positive and negative £ the potential can be approximated as: 

U(Q = U+(t)Q(0 + U-(0Q(-0, (14) 



2n + 1 



£ + 5 (2n+l)/(2n+2A 2n/(2n+1) _ g 2n/(2n+2) 



where 



2ff 

There is a kind of the conservation law for the energy of field £: 



(15) 



1 f T Bz 

^' 2 + U(t)- J dr'—^r')= const, 



(16) 



where n is an arbitrary fixed time moment. The last term appears because U explicitly depends 
on time through z. If dz/dr is positive, which is the case for a contracting body, the value of 
U(£) would in general grow with time. According to assumption made above z linearly grows 
with time as z(r) = 1 + kt, where 

k = {mtcontr^fgY 1 . (17) 
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Figure 1: Oscillations of £i(t), in the case n = 2, k = 0.01, g = 0.01 and initial 
conditions tjq = 1, y' Q = k/2. The amplitude of the oscillations is in very good 
agreement with analytical estimate (|12|) . 



This simple law may be not accurate when t/t con t r > 1, but probably the results obtained are 
not too far from realistic case. 

Equation of motion dSJ) for small oscillations £i can be rewritten as 

£l' + n 2 £i = (is) 

Term £™j n is proportional to k 2 , which is usually assumed small, so in first approximation it 
can be neglected, though an analytic solution for constant O or in the limit of large O can be 
obtained with this term as well. Using expansion ()12p and neglecting a" , we obtain 

1 9 \ 1 L , 3 




Here and in what follows sub-0 means that the corresponding quantity is taken at initial moment 

T = T . 

We impose the following initial conditions 

[y{r = r ) = z(t = T ) = 1 , ^ 
\y'(T = r ) = y' , 

the first of which corresponds to GR solution at the initial moment. The initial value of the 
amplitude derivative ckq can be expressed through y' , which we keep as a free parameter, as: 



{K-y' Q ){g + 2n+lfl 2 . (21) 



To keep £i small the initial value of the derivative y' should be also small. In this case the 
numerical results, shown in figure [lj are in remarkable agreement with eq. f)12|) . For y' = k the 
oscillations seem not to be excited. However, this is an artifact of the approximation used. For 
example, the "source" term in the r.h.s. of eq. (|18p creates oscillations and hence deviations from 
GR with any initial conditions. 

Our primary goal is to determine the amplitude and shape of the oscillations of y. In contrast 
to £ the oscillations of y are strongly unharmonic and for negative and even very small |£| the 
amplitude of y may be very large because y —£/g, according to eq. Q. This feature is well 
demonstrated by the results of numerical calculations shown in figure [2j 

We can estimate the amplitude of the spikes using the energy evolution law (|16|) . Let us use 
it at the moment when the maximum value of £ reaches zero. So at this moment £' = 0. Since 
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Figure 2: "Spikes" in the solutions. The results presented are for n = 2, g = 0.001, 
k = 0.04, and y' Q = k/2. Note the asymmetry of the oscillations of y around y = z 
and their anharmonicity. 



U(0) = 0, the constant in the r.h.s. of eq. (|16|) turns to zero. Now let us go to larger time and 
neglect the oscillating part of £ under the integral. The maximum absolute value of negative £ is 
determined by the equation: 



C / dr' 



-(2n+l) 



1 

2n 



1 



2» 



The value of z(t\) = z\ is found from the condition a = ^ m in 

Z~ {2n+1) -g Zl = \ K -y> \(g + 2n + lf* 

In the limit of small g, g < l/zf n+2 , asymptotically |£ ma z| = 
eq. (|23|) gives: 



z(r) 



i.e. 



2n 



+ - 5 [z 2 (n) - z\t)] . 

(22) 



2n + 1 

2n+2 
1 



Z 



(9/n) 



1/2 



-1 



Vo) 2 (2n + l) 



-l/(3n+l) 



This finally determines 



(ng) 



-1/2 - n 



-1 



This result is much larger than the naive estimate gytl 



2n+2 



1 



1/4 

(23) 

z(Ti)~ n . In the same limit 
(24) 
(25) 

For negative £ potential 



(|15p behaves as C/ ~ ^ 2 /(2g), so the characteristic frequency of the oscillations in the region of 
negative £ is about O ~ i n dimensionless time which corresponds to cj ~ m in physical 

time. Evidently frequency of oscillations of y in this region is the same. It can be shown that the 
calculated amplitude of y (|25p becomes noticeably larger with rising z [16] . 

According to calculations of ref. [14] , harmonic oscillations of curvature with frequency oj and 
amplitude R ma x transfer energy to massless particles with the rate (per unit time and volume): 



QPP^Rl ax u;/(11527r) 



(26) 



The life-time of such oscillations is tr = 48 m 2 pl /uj^. 

In our case the oscillations are far from harmonic and we have to make Fourier expansion 
of the spiky function y(r). To this end we approximate the "spike- like" solution as a sum of 
gaussians of amplitude modulated by slowly varying amplitude B(t), superimposed on smooth 
background A(t): 



R(t) = A(t) + B{t)J2exp 



2a 2 



(27) 
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We assume that a <C t\, that is the spacing between spikes is much larger than their width. The 
Fourier transform of expression (f2"7|) is straightforward but rather tedious (details can be found 
in our work [16]). Finally we find: 

km ~ - 2 E <5 ( w --rJ- ( 28 ) 

1 y \ 

Identifying i? with R m ax = VmaxTo and integrating over frequencies we obtain 

gpp^yla^/^GTTh). (29) 

Time interval ii is approximately equal to 2tt/lj s i ow = 2tt /(Q s i ow my/g), where £l s iow is given by 
eq. (HH). Amplitude y m ax is defined in eq. ([25]) ; parameter k is determined by eq. ([T7j) and 
"coupling constant" g is defined below eq. ([7]). Taking all the factors together we finally obtain: 

n-l (n + l)(7ra+l) 

frp- ^y^ 1 (JO m r») , ( 3 o) 

^ pi^contr \ tcontr J \ f?0 / 

where 

9n-l 7n + l 

C n = (2n + l)2(3n+i) (GAn) 2 ^ 1 ) /(18n). (31) 



It is convenient to present numerical values: g c /mp, ~ 2 • 10 123 and (mtu) 2 ~ 3.6 • 10 93 



where g c ~ 10 _29 g/cm 3 and = m/10 5 GeV. Now assuming that the particle production lasts 
during t ps t CO ntr and taking Qo = Q c , we find the energy flux of cosmic rays produced by oscillating 
curvature: 

qcr « 10~ 24 m^ z n+1 GeVs" 1 cm" 2 . (32) 

This result is a lower limit of the flux of the produced particles. With larger z when the minimum 
of the potential shifts deep into the negative £ region the production probability significantly rises. 
Analytical evaluation in this case is more difficult but simple qualitative arguments indicate to 
it and numerical calculations supports this assertion. 

Recall that our derivation has been done under assumption R > Rq and respectively £o > Qc- 
However if we take qo somewhat larger than q c the results would be approximately correct. 

For m = 10 10 GeV the predicted flux of cosmic rays with energy around 10 19 eV is close or 
even higher than the observed flux. 
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